
install.packages("flexmix")

library(flexmix)


#setwd("~/Dropbox/Job Market Paper/IO Revisions/Replication_Files/Replication_Appendix/A_Fig_10")





data_q<-read.csv(file="data_fig_A10.csv")

years<-unique(data_q$Year)

one<-stepFlexmix(logArea~1| Year,k=1:7,data= data_q)


print("this shows the model w/ 3 is best")

pdf(file="NOT_IN_APPENDIX_BUT_TO_CONFIRM BEST MODEL.pdf")
plot(one)
dev.off()

two<-flexmix(logArea~1|data_q$Year,k=3,data=data_q)

clusters(two)

data_q$cluster<-clusters(two)


unique(subset(data_q,data_q$cluster==1)$Year) #1490-1790
unique(subset(data_q,data_q$cluster==2)$Year) # 1245-1485
unique(subset(data_q,data_q$cluster==3)$Year) # 1100-1240


#############################################




quantile_1Lts<-tapply(data_q$logArea,data_q$Year,quantile,probs=c(.1))

quantile_9Lts<-tapply(data_q$logArea,data_q$Year,quantile,probs=c(.9))

quantile_5Lts<-tapply(data_q$logArea,data_q$Year,quantile,probs=c(.5))





###############################################


pdf(file="Figure_A10.pdf")

plot(years,quantile_1Lts,type="n",ylim=c(3.3,12.3),xlab="Year",ylab="log(Km^2)")


xx1<-c(years[1:27.5],rev(years[1:27.5]))
Y1<-c(quantile_1Lts,quantile_9Lts)[1:27.5]

yy1<-c(quantile_1Lts[1:27.5],rev(quantile_9Lts[1:27.5]))

polygon(xx1,yy1,col="lightgrey")




xx2<-c(years[27.5:77.5],rev(years[27.5:77.5]))

yy2<-c(quantile_1Lts[27.5:77.5],rev(quantile_9Lts[27.5:77.5]))

polygon(xx2,yy2,col="grey")


xx3<-c(years[77.5:139],rev(years[77.5:139]))

yy3<-c(quantile_1Lts[77.5:139],rev(quantile_9Lts[77.5:139]))

polygon(xx3,yy3,col="darkgrey")


lines(years,quantile_5Lts,type="l",ylim=c(3.3,12.3),xlab="Year",ylab="log(Km^2)")
lines(years,tapply(data_q$logArea,data_q$Year,mean))


#lines(years,quantile_1Lts,lty=2)

#lines(years,quantile_2Lts,lty=3)

#lines(years,quantile_3Lts,lty=4)

#lines(years,quantile_4Lts,lty=5)

#lines(years,quantile_6Lts,lty=7)

#lines(years,quantile_7Lts,lty=8)

#lines(years,quantile_8Lts,lty=9)

#lines(years,quantile_9Lts,lty=10)

#lines(years,quantile_5Lts,lty=1)

text(1790,10.5,"90th")
text(1790,3.5,"10th")


dev.off()

